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Abstract. The stability of rotating isotropic spherical stellar systems is investigated by using N-body simulations. Four spher- 
ical models with realistic density profiles are studied: one of them fits the luminosity profile of globular clusters, while the 
remaining three models provide good approximations to the surface brightness of elliptical galaxies. The phase-space distribu- 
tion function f(E) of each one of these non-rotating models satisfies the sufficient condition for stability df /dE < 0. Different 
amounts of rotation are introduced in these models by changing the sign of the z-component of the angular momentum for a 
given fraction of the particles. Numerical simulations show that all these rotating models are stable to both radial and non-radial 
perturbations, irrespective of their degree of rotation. These results suggest that rotating isotropic spherical models with realistic 
density profiles might generally be stable. Furthermore, they show that spherical stellar systems can rotate very rapidly without 
becoming oblate. 
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1. Introduction 

Dynamical instabilities in spherically symmetric stellar sys- 
tems have been investigated for more than four decades. In a 
seminal work, Antonov (1960 1962) used a variational prin- 
ciple to demonstrate that non-rotating spherical models with 
a phase-space distribution function / depending only on the 
energy E are stable to non-radial perturbations if df/dE < 0. 
Subsequent works showed that this condition is also a sufficient 
condition for stability to radial perturbations (Doremus et al. 
119711; Sygnet et al. 119841; Kandrup & Sygnet 119851). In general, 



non-rotating spherical stellar systems are described by distribu- 
tion functions / that depend on both the energy E and the mag- 
nitude of the angular momentum L. In such systems only the 
stability to radial modes can be tested by using the sufficient 
condition df/dE < (Doremus & Feix 1973; Dejonghe & 



Merritt 1988 1. For this reason, numerical simulations have been 



an indispensable tool to investigate the stability of anisotropic 
spherical models. 

Several classes of instabilities have been discovered in non- 
rotating spherical model s with anisotropic velocit y distribu- 
tio ns (e. g., Henon |1973[ M erritt & Aguilar |1985| ; Barnes et 
al. |1986| ; see Merritt |1999| for a recent review). For example, 
models dominated by stars on radial or eccentric orbits can be 
unstable to forming a triaxial bar (Polyachenko 1981; Merritt 
& Aguilar 1985; Meza & Zamorano 1997), while models com- 
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posed mainly of stars on circular orbits can exhibit quadrupole- 
mode oscillations (Barnes et al. 1986). 



In contrast to non-rotating spherical models, comparatively 
little work has been done to investigate the stability of rotat- 



ing spherical stellar systems. Miller & Smith (1980) employed 
several methods to introduce rotation in a spherical n — 3 
polytrope with isotropic velocity distribution. In all cases, they 
found that the addition of rotation does not affect the stability 
of their models; indeed, they found that even rapidly rotating 
syst ems r emain spherically symmetric. More recently, Alimi et 
al. ( 1999| ) showed that rotating isotropic spherical n — 4 poly- 
tropes, which were made to rotate by changing the sign of the 
z-component of the angular momentum for a fraction of the 
particles, were dynamically stable. 



On the other hand, Allen et al. ( 1992 ) reported the existence 
of a "tumbling bar instability" in a set of rotating spherical 
models with different degrees of velocity anisotropy, ranging 
from completely circular to entirely radial (see also Papaloizou 
et al. |1991|; Palmer |l994a[ |l994b[). In particular, they sug- 



gested that the introduction of a small amount of rotation in 
the isotropic spherical n — 2 polytrope induces a bar instability 
in this otherwise stable system. However, Sellwood & Valluri 
(1997) using the same files of initial conditions but a different 
N-body code, showed that these models are actually stable; the 
instability appears do not exist. They suggested that the evolu- 



tion observed by Allen et al. (1992) was, probably, caused by 



an improper treatment of variable time steps in their N-body 
code. 

Most of the previous results have been obtained for rotating 
spherical models with somewhat unrealistic properties. Indeed, 
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most of these works have employed spherical polytropes with 
finite radius (e.g., Alimi et al. 1999) or models with unrealistic 



Table 1. Parameters for the models 



velocity distributions (e.g., Allen et al. 1992). Therefore, they 
cannot be considered as general. To investigate the influence of 
rotation on the stability of spherical stellar systems, it is neces- 
sary to consider models with more realistic density profiles and 
velocity distributions. In this paper, the results of a series of N- 
body simulations for four of such models are presented. These 
isotropic spherical models are: the Plummer (1911) model, the 



Hernquist ( |1990) mod el, the Jaffe ( |1983|) m odel, and the y = 
model (Dehnen 1993; Tremaine et al. 1994 1. Different degrees 



of rotation were introduced in these models by using the so- 



called Lynden-Bell's (1960 1962) demon to reverse the sense 



of rotation along the z-axis of a given fraction of the particles. 
Numerical simulations show that all these rotating models are 
stable, regardless of their degree of rotation. 

This paper is organized as follows. The method employed 
to introduce net rotation in these models and the N-body 
code used to follow their dynamical evolution are described 
in Sect. 2. The main results of the numerical simulations are 
summarized in Sect. 3. Finally, a brief discussion of the results 
is given in Sect. 4. 

2. Numerical simulations 

2.1. Rotating models 

The sp herical models studied in this paper are: the Plummer 
fll91l|) model, 



Ma 1 



An (a 2 + r 2 ) 5 > 2 ' 



the Jaffe fll983| ) model, 
1 Ma 



P = 



An r 2 (a + r) 2 ' 



the Hernquist (199C) model, 
1 



Ma 



P = 



2n r(a + r) 3 ' 

and the y — model (Dehnen 1993 Tremaine et al. 1994), 



aM 



P = 



An (a + r) 4 ' 



(1) 



(2) 



(3) 



(4) 



In all these cases, M is the total mass and a is the scale ra- 
dius. Hereafter, all quantities are expressed in units such that 
M = a = G = 1, where G is the gravitational constant. Some 
parameters for these models are summarized in Table [j], where 
rh is the half-mass radius and t/, is the dynamical time evaluated 
at r h . 

The Plummer model fits the light distribution of globular 



Model 




h 


At 


Plummer 


1.31 


3.32 


0.008 


Hernquist 


2.41 


8.33 


0.02 


Jaffe 


1.00 


2.22 


0.006 


7 = 


3.85 


16.78 


0.04 



phase-space distribution function f(E) can be obtained ana- 
lytically by using th e Edd ington's inversion formula (see e.g., 
Binney & Tremaine 1987 ). For all these non-rotating models, 
the distribution function satisfies the sufficient condition for 
stability df/dE < 0. Therefore, the models are st able to both 
ra dial an d non-radial p erturb ations (Antonov 1962 ; Doremus et 
al. |1971| ; Sygnet et al. |1984| ; Kandrup & Sygnet |l985| ). 

Initial conditions for the simulations were derived from the 
density profile and the distribution function by using the fol- 
lowing scheme. The radial coordinate of each particle is as- 
signed by inverting the equation M(r) = x, where M(r) is the 
total mass inside the radius r and x is a uniform random vari- 
able in the range [0, 1). Then, the coordinates of the position 
vector x are chosen at random from a sphere with radius r. The 
distribution function f(E) and the gravitational potential O(r) 
evaluated at the particle position define the distribution of ve- 
locities, which is sampled by an acceptance-rejection technique 
to assign the modulus of the velocity v of each particle (see 



e.g., Press et al. 1992). Finally, the coordinates of the velocity 
are chosen at random from a sphere with radius v. 

The distribution function f(E) does not depend on the sign 
of the z-component of the angular momentum L z (or equiv- 
alently on the sign of v<j), therefore the models have no net 
streaming. Different amounts of rotation can be introduced in 
these models by reversing the sense of rotation about some 
axis, here taken to be the z-axis , of a given fraction of the parti- 
cles (Lynden-Bell 1960 1962 ). In a non-rotating model, there 
are equal number of particles with going in opposite direc- 
tions, while for a maximally streaming model the sign of 
is the same for all particles. All intermediate cases have vary- 
ing fractions of particles with going in opposite directions. 
This scheme preserves the position and the norm of the veloc- 
ity of each particle, then the systems are put in rotation without 
modifying their total kinetic and potential energy. Therefore, 
the resulting rotating models are also in dynamical equilibrium 
(see Lynden-Bell |1960| ). 

This flipping rule introduces a discontinuity in the distribu- 
tion function across L, = 0, which, in principle, can enhance 



the strength of the bar mode (Kalnajs 1977). However, as it is 
shown below, no signs of bar instabilities are observed in these 
rotating models. Therefore, no special procedure was necessary 
to taper this flipping rule when \L Z \ is small. 

The degree of rotation in these models can be measured by 



clusters (see e.g., Spitzerp^, while the remaining three den- the P™eter (see e.g., Sellwood & Valluri ^997j) 



sity profiles provide good approximations to the surface bright- 



ness of elliptical galaxies (see Dehnen 1993). Several intrinsic 
properties and projected quantities of these models can be ob- 
tained analytically, e.g., the velocity dispersions, the surface 
brightness profile and the mass distribution. In particular, the 



n = 



2f =I I/, 



(5) 



which varies from rj — for a non-rotating model to 77 = 1 for a 
model with all particles orbiting in the same sense around the z 
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Table 2. Parameters for the rotating models. 



0.3- 



1 


C 





0.17 


0.2 


0.20 


0.4 


0.23 


0.6 


0.27 


0.8 


0.30 


1.0 


0.33 



axis. When 77 = 0.5 the system has half the maximum possible 
total angular momentum; in this case, 75% of the particles are 
orbiting in the direct sense and the remaining 25% are retro- 
grade. 

Alternatively, the amount of rotation can also be given in 
terms of the following parameter (see Navarro & White 1993) 



^rot 

K ' 



(6) 



where K is the total kinetic energy and K mt is the rotation ki- 
netic energy defined by 



K n 



1 N 

1=1 1 



mi(Li 



Ltot) 2 



(r; 



■ L tot ) 2 



(7) 



In order to exclude counterrotating particles, the sum in equa- 
tion (Q) is carried out only over particles that satisfy the condi- 
tion (L, • L t ot) > 0, where L tot is a unit vector in the direction of 
the total angular momentum of the system, in this case, along 
the z-axis. Therefore, this parameter measures the kinetic en- 
ergy in the rotational motion around the z-axis. It varies from 
(i=l/6 for a non-rotating model to (i = 1/3, the maximum 
value allowed for the virial theorem, for a model with all parti- 
cles rotating in the same sense around the z-axis. Table || sum- 
marizes the values of rj and the associated values of (i employed 
in the simulations. 

Figure [j] shows the mean streaming velocity for Plummer 
models with different degrees of rotation. To measure this 
quantity, the systems were first projected along the y-axis and, 
then, a slit was placed along the projected x-axis. The slit width 
was set at 77,, the half-mass radius, and the bin size along the 
projected axis was varied such that each bin contains the same 
number of particles (n D in = 1000). In all the rotating models, 
the velocity curve reaches the maximun at R 55 77, and then de- 
creases slowly (models with 77 0.6) or remains practically 
flat (models with < rj 55 0.6). As expected, the mean stream- 
ing velocity of the non-rotating model (77 = 0) is nearly zero; 
the observed fluctuations reflect the discreteness of the models. 
The velocity curves for the other models are very similar. 

2.2. N-body code 

The stability of these models was investigated by using an N- 
body code based on the self-consistent field method described 



by Hernquist & Ostriker ( 1992 1. This scheme consists of solv 




Fig. 1. Mean streaming velocity for Plummer models with dif- 
ferent degrees of rotation 77 as a function of the projected po- 
sition along the x-axis. The systems are viewed such that the 
rotation axis is perpendicular to the line-of-sight. 



gravitational potential using a biortoghonal set of basis func- 
tions. For spherically symmetric systems, it is natural to em- 
ploy spherical harmonics to expand the angular dependence. 
Then, the density and potential expansions become 



p(r) = ^A n i„,p n i(r)Yi, n (6,(p), 

nlm 

<t>(r) = Yj A "im^ni(r)Yi m (e,(p). 



(8) 
(9) 



In practice, these expansions are truncated at some values 7i max 
and / max for the radial and angular functions, respectively. 

The choice of the radial basis functions {p„i, <£„;} is not 
unique; in fact, several sets have been proposed (e.g., Clutton- 



Brock 1973; Allen et al. 1990; Zhao 1996). However, the ef- 



ficiency of this method relies upon the ability to represent the 
density profile of the initial system and its subsequent evolu- 
tion with the first few terms of the basis set. Therefore, it is 
generally desirable to select a set of basis functions such that 
their lowest order terms provide a good approximation to the 
density profile of the system under study. 



In these simulations, the Clutton-Brock (1973) basis set 



is used for the Plummer model, while the Hernquist-Ostriker 



(1992) basis set is employed for the other three models. The 
zeroth-order terms of these basis sets are, respectively, the 
Plummer model and the Hernquist model. The y = model can 
be recovered by a linear combination of only two terms of the 
Hernq uist-Ostriker basis functions (see Hernquist & Ostriker 



ing the Poisson equation by expanding the density and the 



1992J). Therefore, these basis sets provide exact representations 
for three of the models here studied. On the other hand, the 
Jaffe model can only be approximated by using a finite number 
of terms of the Hernquist-Ostriker basis set; in this case, good 
accuracy c an be obtained with ~ 10 terms (see e.g., Meza & 
Zamorano 1997 ). 

All simulations employed N = 10 5 equal-mass particles 
and potential expansions up to n max = 10 for the radial func- 
tions and Z max = 2 for the angular functions. Nonzero m 
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Fig. 2. Evolution of the axis ratios c/a and b/a, with a > b > c, for the y — models with different degrees of rotation. The 
axis ratios were measured at radius r m = 8, which encloses ~ 70% of the total mass. Time has been normalized to the half-mass 
dynamical time, f/i = 16.8. 



terms were also included to allow the development of any non- 
axisymmetric instabilities. Time integration was performed us- 
ing a second order scheme with a fixed time step Af, given by 

1 9 

x, + i = x, + Afv, + -Af 2 a ; , (10) 

v i+ i = v,- + i At (ai + a /+ i), (11) 

wher e the subscript identifies the iteration (see e.g., Hut et al. 
1995). The total elapsed time for all these simulations was 
Tend = 50th, where f/, is the dynamical time evaluated at the 
half-mass radius of the respective model. The time step used for 
each model is shown in Table [j]. With these parameters, the to- 
tal energy was conserved to better than 0.01 % for the Plummer, 
Hernquist and y = models, while for the Jaffe models the 
conservation of energy was only about 4%. 

With the adopted value Z max = 2 for the potential ex- 
pansions, these simulations are designed for searching radial 
(m = 0), lopsided (m =1), and bar (m = 2) instabilities. A test 
particularly sensitive to the bar instability is the evolution of 
the axis ratios of the particle distribution. For each simulation, 
the axial ratios inside a given radius were obtained by using an 
iterative algorithm similar to that used by Dubinski & Carlberg 
( 1991[ ). In this scheme, initial values for the inertia tensor 

< 12 > 



are computed for all particles inside a sphere of radius r m . The 
eigenvalues and eigenvectors of 7 (J provide an approximation to 



the axis ratios and the orientation of the fitting ellipsoid. Then, 
the modified inertia tensor is evaluated only for particles in- 
side that ellipsoid, which gives an improved approximation to 
their axis ratios and orientation. This process is repeated un- 
til the axis ratios converge to a value within a pre-established 
tolerance. 

3. Results 

A total of 24 simulations were performed. The main result of 
these simulations is that all the rotating models are dynamically 
stable. This is illustrated by displaying the evolution of some of 
them. The results for other models are similar. 

The evolution of the axis ratios b/a and c/a, where a > b > 
c, for the y = models with different degrees of rotation is 
shown in Figure |[ These axis ratios were measured at radius 
r m = 8, which encloses ~ 70% of the total mass. For all values 
of 77, the axis ratios remains essentially equal to their initial 
values. The fluctuations reflect the discreteness of the models 
and are consistent with the number of particles employed in 
the simulations. Similar behavior is observed for the axis ratios 
measured at other radii. 

A more sensitive test for the existence of possible insta- 
bilities is provided by the analysis of the individual expansion 
coefficients A„;,„, which is straightforward in the self-consistent 
field method because the coefficients are evaluated at each time 
step in order to compute the forces. The evolution of the am- 
plitude of several expansion coefficients A„\ m for the rotating 
Hernquist model with r\ — 1 is shown in Figure |J These 
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Fig. 3. Evolution of the amplitude of several expansion coefficients A„i,„ for the Hernquist model with rj = 1 . The coefficients 
A„i m were computed using the Hernquist-Ostriker basis set, where n refers to the radial functions, and I and m to the angular 
functions (spherical harmonics). The logarithm of the absolute value of the coefficients is shown. Time has been normalized to 
the half-mass dynamical time, ?/, = 8.33. 



coeffic ients were calculated using the Hernquist & Ostriker 
( 1992 ) basis functions. As expected, the zeroth-order coeffi- 
cient Aooo - 1, while the higher order coefficients are closer 
to zero and fluctuate around their initial values with amplitude 
consistent with the noise induced by the number of particles 
used in these simulations. In particular, there are no signs of 
radial (m = 0), bar (m = 2), or lopsided (m =1) instabilities. 
Other coefficients show the same qualitative behavior. 

4. Conclusions 

The stability of four rotating isotropic spherical models was in- 
vestigated by using N-body simulations. The density profiles 
of these models provide good approximations to the mass dis- 
tribution of globular clusters and elliptical galaxies. Different 
degrees of rotation were introduced in these models by revers- 
ing the sense of rotation along the z-axis of a given fraction 



of the particles (Lynden-Bell 1962). Simulations show that all 



these rotating models are dynamically stable, irrespective of 
their degree of rotation. No signs of radial, lopsided, or bar in- 
stabilities were observed. If some of them exist their growth 
rate is larger than 50th, the time elapsed for these simulations. 

These simulations show that spherical stellar systems can 
rotate very rapidly without becoming oblate. This result con- 



trasts with the suggestion of Alimi et al. ( 1999 ) that there do 
not exist spherical stellar systems in fast rotation (i.e., with 
/j. Z 0.1). However, their conclusions were based on a series of 
simulations for spherical n = 4 polytropes with isotropic veloc- 
ity distributions, which were made to rotate by using a proce- 
dure that modifies their initial velocity anisotropy. Therefore, 
their conclusions are not strictly applicable to the isotropic 
models studied in this paper. But, they could still be valid for 
anisotropic spherical models. In fact, there are several other 
ways to construct rotating spherical models. For example, one 
could construct a system containing only circular orbits and 
then reverse a half of them (see Lynden-Bell 1960). Clearly, 
such a model would exhibit stronger rotation that the models 
here analyzed and surely might be unstable. 
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